First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other 
molecules of interest for organic photovoltaic applications 
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We evaluate the performances of ab initio GW calculations for the ionization energies and HOMO- 
LUMO gaps of thirteen gas phase molecules of interest for organic electronic and photovoltaic ap- 
plications, including the Ceo fullerene, pentacene, free-base porphyrins and phtalocyanine, PTCDA, 
and standard monomers such as thiophene, fluorene, benzothiazole or thiadiazole. Standard Go Wo 
calculations, that is starting from eigenstates obtained with local or semilocal functionals, signifi- 
cantly improve the ionization energy and band gap as compared to density functional theory Kohn- 
Sham results, but the calculated quasiparticle values remain too small as a result of overscreen- 
ing. Starting from Hartree-Fock-like eigenvalues provides much better results and is equivalent to 
performing self-consistency on the eigenvalues, with a resulting accuracy of 2-4% as compared to 
experiment. Our calculations are based on an efficient gaussian-basis implementation of GW with 
explicit treatment of the dynamical screening through contour deformation techniques. 

PACS numbers: 71.15.-m,71.15.Ap,71.15.Qc,71.20.Rv 
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I. INTRODUCTION 

The flexibility in the synthesis of novel molecules and 
polymers is an important advantage of organic photo- 
voltaics as compared to the inorganic routei*^. Despite a 
rather limited quantum efficiency, the possibility to tai- 
lor the solubility, cristallinity and electronic properties 
of the building molecular units is offering much means 
to improve on the actual best cells, such as those based 
on the combination of acceptor fullerene derivatives and 
derivatives of polythiophene as donors^. In particu- 
lar, it has been shown that there are strong correlations 
between the "band offsets" at the donor/acceptor inter- 
face and the open circuit voltage or the driving force 
for separating the hole and electron of the photoinduced 
excitons^. The ability to tune the electronic affinity and 
ionization energy of the donor and acceptor molecules, 
under the constraint that sun light absorption should be 
kept as large as possible, is a current and intense field of 
research^"— . There is therefore much interest in devel- 
oping efficient quantum simulation methods allowing to 
provide the spectroscopic and optical properties of stan- 
dard molecules with both a reasonnable computer cost 
and accuracy. 

For isolated molecules, an excellent trade-off between 
computer cost and accuracy for the calculations of the 
ionization energy and electronic affinity can be found 
with the so-called ASCF approach using hybrid function- 
als such as PBEO and B3LYP obtained by admixture of 
a fixed amount of Fock exchang o^^'^^ . However, these 
techniques cannot be used for extended systems such as 
bulk semiconductors, molecules deposited on a surface or 
in solution, and the percentage of Fock exchange needed 
for obtaining good results with these functionals is ex- 
pected to change from isolated molecules to bulk sys- 



tems. For the same reasons, the "Kohn-Sham" ioniza- 
tion energies, electronic afflnities and band gaps as ob- 
tained from the eigenvalues of the Hamiltonian may be 
certainly improved with hybrid functionals as compared 
to (semi)local ones, but again the amount of Fock ex- 
change needed to get accurate results may change from 
one system to another. 

A technique based on many-body perturbation theory 
(MBPT), namely the GW approximationi^"— , has shown 
excellent results for the evaluation of the band edges 
and band gaps of extended bulk systemsii. Distinct 
from the perturbative techniques developed by the quan- 
tum chemistry community to build up correlations from 
the Hartree-Fock solutionis, such an approach is gen- 
erally derived from functional derivative technique a^*^'^^ 
yielding an exact (non-perturbative) set of self-consistent 
(closed) relations between the one-body Green's func- 
tion G, the polarizability P, the dynamically screened 
Coulomb potential W, the "exchange and correlation" 
self-energy S and the so-called vertex corrections P, 
which is related to the variation of the self-energy with 
respect to an external pertubation. In practice, neglect 
of vertex corrections leads to the so-called "GW" ap- 
proximation for the self-energy which can be loosely de- 
scribed as a generalization of the Hartree-Fock method 
by replacing the bare Coulombian with a dynamically 
screened Coulomb interaction. The ingredients needed 
to proceed through the GW calculations pave further the 
way to Bethe-Salpeter calculations^^ aiming at explor- 
ing optical absorption spectra as an alternative to time- 
dependent DFT. While decades of expertise exist for ap- 
praising the performances of the GW approximation in 
the case of extended bulk systems, the application of such 
MBPT approaches to organic molecules in the gas phase, 
and in particular molecules of interest for photovoltaic 
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applications^Sr— , remain extremely scarce, a situation 
that can be mostly attributed to the associated compu- 
tational cost for molecules such as fullerene derivatives or 
porphyrins containing several dozens of atoms. As a re- 
sult, an understanding of the merits of such an approach 
in the case of organic molecular systems, as compared to 
well-established quantum chemistry techniques, is still in 
its infancy. 

We present in this work a GW study of the quasi- 
particle properties of thirteen of the most standard 
molecules involved in organic electronic and photo- 
voltaic devices, including the Cgo fullerene, the free- 
base 21H,23H-porphine (H2P), tetraphenylporphyrin 
(H2TPP), and phtalocyanine (HzPc), and the 3,4,9,10- 
perylene tetracarboxylic acid dianydride (PTCDA) (see 
Fig. [1]). We also study anthracene, tetracene, and pen- 
tacene, 7r-conjugated molecules of interest for organic 
electronics, even though not as such for optical appli- 
cations, and for which experimental band gap data arc 
available. Finally, the tiophene, fluorene, benzothiazole, 
2,1,3-benzothiadiazole and 1,2,5-thiadiazole monomers, 
building blocks of common donor polymers, are also 
invcstigated^^^. Our results suggest that while the 
standard non-self-consistent GqWo calculations based on 
Kohn-Sham eigenstates with (semi)local functionals cer- 
tainly improves on the DFT results, the GqWq ioniza- 
tion energy and HOMO-LUMO gap remain underesti- 
mated as compared to experiment. A simple partial self- 
consistency on the eigenvalues only, or the use of Hartree- 
Fock-like eigenvalues in a one-shot GqWo calculation, al- 
lows to obtain much improved results. We show in par- 
ticular that these simple schemes lead to an average error 
of ~0.3 eV for the ionization energies and 0.1-0.2 cV for 
the band gaps. 

Our paper is organized as follows. In section (II), we 
briefly describe our implementation of the GW formalism 
within a gaussian-basis, including details about the eval- 
uation of the Coulomb matrix elements. In section (III), 
our results for the ionization energy and HOMO-LUMO 
gap of selected molecules are presented and compared to 
existing experimental results. The importance of a sim- 
ple self-consistency on the eigenvalues is discussed. Sec- 
tion (IV) describes a simplified non-self-consistent ap- 
proach based on an approximate perturbative Hartree- 
Fock starting point for building the Green's function and 
screened Coulomb potential. We conclude in section (V). 



II. METHODOLOGY 

Our code is based on a gaussian-basis implementation 
of the GW formalism and builds on a previous implemen- 
tation of calculating the inverse dielectric matrix using 
numerical strictly localized orbitals^l. To avoid dealing 
with numerical basis, the present implementation now ex- 
pands the needed two-point operators (bare and screened 
Coulomb potentials, susceptibilities, etc.) on an "auxil- 
iary" gaussian basis composed of one-center atomic-like 



orbitals, with real spherical harmonics for the angular 
part and a radial dependence composed of gaussian func- 
tions. The use of such an auxiliary basis, commonly im- 
plemented in several DFT quantum chemistry codes to 
express the charge density for ground-state or excited- 
state^ calculations, allows to greatly speed up the eval- 
uation of e.g. the Coulomb matrix elements. We discuss 
these points in the following subsections. 



A. General formalism 

With the notations of Ref. [2^ we introduce for any 
two-point function /(r,r') the < / > and [/] matrices in 
the auxiliary basis related through: 

[fUu^ J I dvdv' A.*(r)/(r,r')Kr') 

where fi and v are elements of the basis and S is the 
overlap matrix. The standard Dyson equation relating 
the dynamically screened Coulomb potential W{u!) to the 
bare Coulomb one (v) can then be written: 

< W{uj) >^<v> + <v> [x°(t^)] < W{uj) > 
with x° the unscreened (free-electron) susceptibility: 



[x°(w)]^^. = H < '^^i^i'^j >< '^ji^i'^' > 

spins i j 

X ( 

\UJ + Ei ~ Ej + 10 Ld — Ei + Ej — 10 J 

where S ~ 0^ . The input {(f>i,Ei) are one-body eigen- 
states and related eigenvalues traditionaly taken as the 
Kohn-Sham solutions of a ground-state DFT calculation. 
In the present paper, we start with a standard DFT/LDA 
calculation but as discussed below, this may not consti- 
tute the best starting point for molecular systems. The 
knowledge of the dynamical screened Coulomb potential 
W{ijj) allows to build the non-local and energy dependent 
self-energy operator S, which accounts for exchange and 
correlation in the present quasiparticle formalisnJ^ and 
reads: 

S°^(r,r'|E) = 7^ / (iwe"^'"G(r,r'|£; + cj)iy(r,rV) 
27r J 

G{r, rV) = J2 <Pn{vWn{^')l{u, - e„ ± i5) 

n 




FIG. 1. (Color online) Symbolic representation of (a) 21H,23H-porphine (H2P), (b) tetraphenylporphyrin (H2TPP), (c) phtalo- 
cyanine (H2PC) (d) 3,4,9, 10-perylene tetracarboxylic acid dianydride (PTCDA), (e) thiopliene, (f) fluorene, (g) benzothiazole, 
(h) 2,1,3-benzotliiadiazole and (i) 1,2,5-thiadiazole. Small white atoms are hydrogen atoms, grey atoms are carbon atoms while 
red/blue/yellow atoms are oxygen/nitrogen/sulfur atoms respectively. 



where the time-ordered Green's function G is again built 
from the ((j)i,Si) eigenstates. The sign of the 6 infinitesi- 
mal insures that the occupied (unoccupied) states corre- 
spond to poles in the fourth (second) quadrants. Again, 
the choice of the "best" input {(j)i,ei) for the building of 
G will be discussed below. 

This implementation is formally equivalent to that 
of Ref. except that we go beyond the plasmon-pole 
model and proceed with the explicit calculation of the fre- 
quency integral for the correlation part of the self-energy, 
^GW ^ j^Giy _ j^^^ ^j^-j^ ^i-^^ Pq^Jj operator. We use 

contour deformation techniques with an integration along 
the imaginary axis complemented by the evaluation of the 
poles in the first and third quadrant for states away from 
the band edgeai^i^: 



The first order perturbation theory self-energy cor- 
rection to the DFT Kohn-Sham eigenvalues is extrapo- 
lated to the quasiparticle energies by a Taylor expansion, 
namely: 



£„ + Z„ < 0„|I]'^'^(£„) - > 



LDA\ 



where Zr, is the renormalization factor defined as: 



i/z„ = i- [as^^/ai?]^^^^. 

with (e„, (f>n) the LDA Kohn-Sham eigenvalues and eigen- 
states in the present case. 



i:^'^{ry\E)=Y,Ur)Kir')Vn{r,r'\E) 



with, introducing W = W ^ v, Ep the Fermi level, and 
9 the Heaviside step function: 



V„(r,r'|S) = W{r, r'|e„ - E) [e{E - £„) - OiEp - e„)] 
'■+°° duj E-En 
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A change of variable allows to fold the smooth func- 
tion W{iuj) onto the finite [0, 1] interval where Gaussian 
quadrature with as little as 12 gaussian points is suffi- 
cient to reach convergency. An analytically integrable 
tail is added/subtracted to avoid instabilities with the 
integrand for w — >■ when i? — ^ e„. 



B. Gaussian basis 

The auxiliary basis used to expand the two-point func- 
tions reads: fi{r) = exp(— ar^)r'i?"'(f) with R\"'{f) the 
real-spherical harmonics and (f) the angular components 
of the r-vector. It is computationally more efficient to 
work with the _R™ instead of the standard YJ™ complex 
harmonics with the following relation: 



[Yr""{f) - (-l)"yr(^)] /V2 (m < 0) 



'{f)] /V2 (m > 0) 
(m = 0) 



The products r^R]^{f) yield the standard expressions 
{x,y,z,xy,yz, i^-?/^, etc.) for the p,d, etc. orbitals (within 
constant factors). We briefly recall that the main ad- 
vantage of a gaussian radial part (as compared to nu- 
merical or Slater-type orbitals) is that the product of 
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two gaussians centered on atoms A and B with de- 
cay coefficients ai and a2 yields a gaussian centered 
on C = {ctiA + a2B)/{al + 02) with a decay constant 
7 = aia2/{cti + ct2)- Further, the r^R!p{f) can easily 
be "shifted" from one center to another with for sake of 
illustration; 

(x - xA){y - va) xc){y - yc) 

+ {yc - yA){x - xc) 
+ {xc - XA){y - yc) + constant, 

showing that a dxy orbital centered on A can be easily 
expressed as a function of {s,p) and dxy orbitals centered 
on C . Such trivial expressions allow to express multi- 
center overlaps in terms of one-center integrals. 

In the present work, our calculations start with a DFT 
calculation of the structural and electronic properties of 
the molecules of interest using the Siesta package^. We 
use a "doublc-C+polarization" (DZP) basis^ and stan- 
dard norm-conserving pseudopotentials. Since the Siesta 
package uses "numerical" orbitals, we first fit the numer- 
ical radial part by up to five contracted gaussians^ in 
order to exploit the relations briefly sketched above. As 
such, both the "ground-state" DFT basis and the auxil- 
iary basis are based on gaussians. Beyond the analycity 
of the gaussian basis, our choice was also motivated by 
the possibility of using eigenstates generated by standard 
chemistry codes with all electron approaches and/or hy- 
brid functionals, providing for some systems possibly a 
better starting point for MBPT calculations (see discus- 
sion below) . We labeled our code "Fiesta" as an attempt 
to extend the "Siesta" package to excited state proper- 
ties. 

Contrary to the planewave case, the auxiliary ba- 
sis for the two-point response functions is larger than 
the ground-state basis. Following Kaczmarski and 
coworkers^, we typically adopt for first raw elements 
such as carbon, nitrogen and oxygen, 4 s,p,d sets of gaus- 
sian orbitals, that is 36 orbitals per atom, while 3 s,p,d 
gaussian orbitals are sufficient for hydrogen. We show 
below that such a basis is large enough for the stud- 
ied organic systems. In the case of sulfur, /-channel or- 
bitals are added. With such a basis, a typical GqWo 
calculation with full dynamics for our largest molecule 
(H2TPP) can be performed within one day on a sin- 
gle standard processor. Better timings and scaling may 
be obtained upon implementing the recently introduced 
techniques allowing to avoid summation over the conduc- 
tion state a^'^''^^" — , or techniques decoupling the sum over 
valence and conduction states^, even though the num- 
ber of unoccupied states is rather limited with standard 
DZP or larger TZDP basis. 

The choice of the "optimal" a-coefficients, controlling 
the localization of the basis orbitals, is a difficult ques- 
tion. Auxiliary basis have been implemented in many 
quantum chemistry codes in order to fit the charge den- 
sity and speed up the calculation of the Coulomb inte- 
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3 in 0.2 3.2 


6.83 
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4 in 0.15 ^ 3.2 


6.89 


6.15 


6.52 


4.74 


7.40 


4.47 


5 in 0.15 3.2 


6.82 


6.06 


6.56 


4.77 


7.29 


4.36 


5 in 0.15 3.5 


6.83 


6.08 


6.51 


4.75 


7.28 


4.33 



TABLE I. Evolution of the ionization (IE) and band gap ener- 
gies of selected molecules as a function of the carbon auxiliary 
basis, changing the number (ng) of gaussians per 1-channel, 
the amin and amax coefficients. Results are in eVs. 



grals. The coefficients of the charge density on the auxil- 
iary basis are optimized using " identity rules"— but not 
the decay coefficients in the exponentials. Years of exper- 
tise in the quantum chemistry community yielded reliable 
auxiliary basis for the periodic table and numerous tests 
have shown that high precision can be obtained with such 
basis provided that they be sufficiently large. 

Since the auxiliary basis must project onto products 
of Kohn-Sham orbitals, optimized basis for all-electron 
calculations cannot be straightforwardly used for GW 
calculations starting from ground-state calculations with 
pseudopotentials. The same guiding lines can however 
be followed. We adopt in particular the idea of a "tem- 
pered" basis^"— suggesting that it is better to generate 
a chain of a parameters such that: aj+i/a, = constant, 
rather than spreading them uniformly between amin and 
Cimax- Such & schemc hinges on the facts that the over- 
lap of two gaussian orbitals is a function of their alpha 
coefficient ratio and that maintaining a constant overlap 
between "adjacent" gaussians allows to better span the 
associated Hilbert spaced. As such, the a„iin , ctmax and 
number of gaussian per 1-channel being chosen, the other 
gaussian coefficients are automatically generated. 

We adopt the basis proposed by Kaczmarski and 
coworkers^, that is namely gaussians with localization 
parameters of (0.2,0.5,1.25,3.2) a.u. for the {s,p,d) 
channels of C, O, and N atoms, and gaussians with 
q:=:(0.1,0.4,1.5) a.u. for hydrogen. As shown in Table I in 
the case of anthracene, II2P porphyrin and Ceo, changing 
the amin and amax values, or increasing the number of 
gaussians in the basis, does not change significantly the 
results. The case of Cgo shows however that reducing the 
number of gaussians to 3 per 1-channcl yields a significant 
error on the band gap. We will show below that the re- 
sults obtained with the present implementation compares 
rather well with previous available calculations based on 
another gaussian basis, planewaves (PWs) or combina- 
tion of gaussians, PWs and Wannier functions. 

We conclude this section related to the auxiliary ba- 
sis by mentioning an important numerical aspect related 
to the overcompletness of the generated non-orthogonal 
gaussian basis. While the basis on a given atom can be 
easily orthogonalizcd using e.g. a Gram-Schmidt proce- 
dure, the overlap between the most diffuse orbitals on ad- 



5 



jaccnt atoms tend to be also rather large yielding an over- 
lap S matrix "nearly singular". Following the strategy 
developed in the case of product-basi; 



Ionization energy 



,39.44 



we rotate our 



basis to that of the eigenvectors of the overlap S-matrix 
from which we remove the eigenvectors with eigenvalue 
smaller than typically 10^^. In the present case of auxil- 
iary basis, such a truncation docs not reduce significantly 
the size of the basis, but avoid the potential numerical 
instability associated with inverting the nearly-singular 
S-matrix and the amplification of errors associated with 
the < V >= S'^^[u]S'^^ transformation (see above). The 
cost of rotating the Coulomb and < (/)i|/3|(/)j > matrix 
elements from the original one-center auxiliary basis 
to the (filtered) S-eigenvectors basis scales as and 
represents a marginal part of the CPU time. 



C. Coulomb matrix elements 

An important ingredient is the evaluation of the 
Coulomb matrix elements between two auxiliary basis 
orbitals localized on two different atoms. Exploiting the 
properties of the Fourier transform (FT) of gaussian- 
based orbitals, namely: 



FT 



(f) = Ce-^''g'i?r(g) 



(1) 



with 7 = l/4a and (f , q) the angular components of the 
(r,q)-vectors in direct and reciprocal space respectively 
(C is a constant), the Coulomb matrix elements reduce to 
a sum of terms built from the product of one-center over- 
laps of three real-spherical harmonics < R^R^ \R\^ > 
(related to Gaunt coefhcicnts with |Z — Z'| < L < {I + 1')) 
times radial integrals 1(1, r;L) of the form: 
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G„W„(LDA) 
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GoWo(HF) 
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GW 




fit on GW 



6.5 7 7.5 8 8.5 9 
Experimental ionization energy [eV] 

FIG. 2. (Color online) Experimental and theoretical ioniza- 
tion energies in electronvolts. Red circles: experimental val- 
ues; light blue triangles up: LDA Kohn-Sham HOMO energy; 
green squares: non-self-consistent Go Wo (LDA) value; black 
diamonds: GW value with self-consistency on the eigenval- 
ues; green stars: non self-consistent GoWo(HFdia9) (see text). 
The black dashed line is a least-square fit of the GW results. 
The figure has been formatted so as to preserve the same 
physical scale on both axis. 



I{IJ';L)^ / dge-«9^g^J,(-/3g2) 
Jo 

The < > factors are pretabulated. The os- 

cillatory behavior of the Bessel function of the first kind 
Ji, makes the direct numerical evaluation rather unsta- 
ble. We prefer to notice that I{1, 1'; L) is straighforwardly 
related to the iFi confluent hypergcometric functionSi^ 
which, for the needed {1,1') values, can be expressed in 
terms of simple hmctions such as the error function (erf) 
or the Dawson integral: F(z) = exp(-z^) crfi(z)/2, 
with erfi(z) = erf(iz)/i, for which rapidly convergent se- 
rial expressions exist^. This is an important advantage 
of the auxiliary basis approach that the evaluation of the 
off-site Coulomb matrix elements is not a costly part of 
the present GW implementation. 



III. RESULTS 

A. Ionization energies 

We start by exploring the ionization energy of our se- 
lected molecules. While experimental data for the elec- 
tronic affinity of molecules in the gas phase are scarce, 
accurate measured ionization energies are much more 
common^. Experimental ionization energies are repre- 
sented by red circles in Fig. [2] and are given in the last 
column of Table II. The DFT-LDA ionization energies, as 
obtained from the opposite sign of the Kohn-Sham high- 
est occupied (HOMO) energy level, are clearly much too 
small, with an average error of f .83 eV or 23% (see blue 
triangles in Fig.[2|). Very similar results are obtained us- 
ing the HOMO energy value as obtained with a scmilocal 



6 



Ionization energy 




LDA-KS 


GoWo(LDA) 


GW 


GoWo(HFdiag) 


Experiment 


anthracene 


5.47 


6.89 


7.06 


7.03 


7.4" 


tetracene 


5.15 


6.37 


6.51 


6.48 


6.97" 


pentacene 


4.94 


5.98 


6.12 


6.08 


6.6" 


Ceo 


6.37 


7.28 


7.41 


7.41 


7.6" 


PTCDA 


6.65 


7.57 


7.68 


7.67 


8.2'' 


HaP 


5.64 


6.55 


6.70 


6.72 


6.9" 


HaTPP 


5.40 


6.09 


6.20 


6.24 


6.4" 


HaPc 


5.56 


6.08 


6.10 


5.93 


6.4^= 


thiophene 


6.15 


8.37 


8.63 


8.64 


8.8" 


fluorene 


5.92 


7.44 


7.64 


7.64 


7.9" 


benzothiazole 


6.33 


8.20 


8.48 


8.50 


8.8" 


thiadiazole 


7.22 


9.65 


9.89 


9.90 


10. 1'' 


benzothiadiazole 


6.55 


8.31 


8.56 


8.57 


9.0" 


MAE 


1.83 (23%) 


0.47 (6%) 


0.30 (3.8%) 


0.31 (4.0%) 





TABLE II. Ionization energies in eV as obtained from the Kohn-Sham eigenvalues (LDA-KS), from non-self-consistent 
GoWo(LDA) calculations, from a GW calculation with self-consistency on the eigenvalues (GW), and from a non-self-consistent 
Go Wo calculation starting from Hartree-Fock-like eigenvalues (GoWo(HFdiag), see text). MAE is the av erag e mean error in eV. 
The average error in percent as compared to the experiment is indicated in parenthesis. "Ref. |47|. ''Ref. [201 '^Ref. |4^. ''Ref. |49|. 



functional such as PBE^. 

We now turn to GoWo(LDA) calculations, that is non- 
sclf-consistent calculations with the Green's function and 
screened Coulomb potential directly built from the LDA 
Kohn-Sham eigenstates and eigenvalues. The analysis of 
the results (column 3 Table I and green squares in Fig. (2]) 
shows that the ionization energies are greatly improved, 
with an average error of 0.47 eV, that is a much reduced 
6% error. 

Even though in much better agreement with experi- 
ment than LDA or PBE, the discrepancies are still size- 
able. As shown below, part of the problem originates 
in that the "starting" LDA HOMO-LUMO gap is dra- 
matically too small for isolated molecules, inducing a 
large overscreening. To avoid using some arbitrary scis- 
sor operator to open the HOMO-LUMO gap in calculat- 
ing the susceptibility, we rather perform a restricted self- 
consistency by reinjecting the corrected eigenvalues in G 
and W up to convergency. As a matter of fact, no more 
than three or four iterations are needed to reach conver- 
gency within 0.01 eV. Such an approximation is labeled 
GW in the following. This is not a full self-consistent ap- 
proach as the eigenstates are not updated, with the ad- 
vantage that the computational cost keeps reasonnable. 
Full self-consistency without vertex corrections is still de- 
bated and seems to yield for small molecular systems re- 
sults that are not as good as GoWq non self-consistent 
runs^. 

The analysis of the results (fourth column Table II and 
black diamonds in Fig. [2|) clearly shows that the self- 
consistency on the eigenvalues improves the results for 
the ionization energy, reducing the average error from 
0.47 eV (6%) to 0.30 eV (or 3.8%). Such a discrepancy 
is still sizeable but much better than the one obtained 
from the LDA Kohn-Sham HOMO energy. An interest- 
ing observation is that the final GW ionization energies 



gather much closer to a straight line (dotted black line on 
Fig. parallel to the first diagonal (red "experimental" 
line) than the LDA data which are much more spread. 
On a pragmatical point of view, this means that the band 
offset between two molecules will strongly benefit from 
cancellation of errors in GW as compared to LDA. In 
particular, the remaining error ('^ 0.2 eV) on the ion- 
ization energy for Cgo, the most standard acceptor, is 
nearly identical to the error on the ionization energies 
of porphyrins and phtalocyanines, which are commonly 
used donors. We now show that self-consistency, even 
though limited to updating the eigenvalues only, has an 
even larger effect on the magnitude of the HOMO-LUMO 
gaps. 



B. HOMO-LUMO gaps 

Due to the lack of experimental values for the elec- 
tronic affinity, experimental quasiparticle HOMO-LUMO 
gaps (red circles in Fig. [3|) are scarce so that we plot our 
results as a function of our "best" calculated HOMO- 
LUMO gaps, namely the GW ones. In the case of Cgo, 
anthracene, tetracene, and pentacene for which exper- 
imental data are available, we observe as expected that 
the LDA HOMO-LUMO gap (blue triangles) is too smaU. 
This is well known in the case of bulk semiconductors but 
here the discrepancy is much larger, with an average error 
of - 4.1 eV or 71%. 

The Go Wo (LDA) HOMO-LUMO gaps (green squares) 
significantly improves with respect to LDA. Compar- 
ing to available Go Wo (LDA) data for this class of aro- 
matic molecules, our calculated 6.15 eV HOMO-LUMO 
gap for anthracene compares well with the 5.97 eV val- 
ues of Niehaus and coworkers, despite the differences 
in basis and the treatment of dynamical effects^. Our 
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GoWo(LDA) 4.79 cV and 4.23 eV HOMO-LUMO gaps 
for the H2P and H2TPP free-base porphyrins respec- 
tively compare further well with the 5 eV and 4.39 eV 
planewave results of Palummo and coworkers^^. Simi- 
larly, our GoWo(LDA) 4.44 eV band gap for Ceo is in 
good agreement with the real-space grid formulation of 
Tiago and coworkers^! yielding a band gap of 4.36 cV. 
Such comparisons certainly underline the reliability of 
the present gaussian-basis implementation. Our 4.53 
eV band gap for PTCDA is however smaller than the 

4.9 eV band gap found with a previous planewave GW 
calculatioit22i56^_ 

Overall, we remark a systematic underestimation of 
the GoWo(LDA) HOMO-LUMO gap with respect to the 
experiment, with an average error for our test molecules 
of ~0.75 eV or 13%. This contrasts with the case of bulk 
systems for which the results of GoWo(LDA) arc gener- 
ally in much better agreement with experimental values. 
Such a behavior can be analyzed by noticing that build- 
ing the polarizabilities and screened Coulomb potential 
with LDA eigenvalues, that is in particular with dramati- 
cally too small HOMO-LUMO gaps, leads to a significant 
ovcrscrecning. This induces too large a correlation cor- 
rection "G(W-V'^)" to the Hartree-Fock HOMO-LUMO 
gap, that is too small a HOMO-LUMO gap. 

Even though much better than the Kohn-Sham 
HOMO-LUMO gap obtained with e.g. the B3LYP 
fimctional^ (see empty down triangles in Fig. 3), it is 
desirable to improve the results. Following the simple 
scheme introduced above, performing self-consistency on 
the eigenvalues in G and W, the GW HOMO-LUMO gap 
is further increased to reach much better agreement with 
experiment. The MAE is now reduced to 0.22 eV or 3.8% 
for our four test molecules. In the case of Ceo, which is 
the most standard acceptor in organic photovoltaic cells, 
the excellent agreement with experiment for the band 
gap value is rather satisfactory. It is interesting to note 
further that the MAE of 0.22 eV for HOMO-LUMO gaps 
is close to the 0.3 eV MAE obtained for the ionization 
energies, suggesting that the electronic affinity is quite 
well reproduced on the average. 



IV. A SIMPLE NON-SELF-CONSISTENT G„Wo 
APPROACH BASED ON HARTREE-FOCK-LIKE 
EIGENVALUES. 

We conclude this study by exploring a simple non- 
self-consistent Go Wo scheme starting from an "ansatz" 
Hartree-Fock (HF) calculation obtained by removing the 
exchange-correlation contribution to the LDA eigenval- 
ues and adding the diagonal part of the exchange opera- 
tor in the LDA basis, namely: 



HOMO-LUMO gap 



."HF" 



LDA 



< 



kLDA\ 



LDA\iLDA 



10. 



where S^, and Vj^'J^^ are the Fock and (semi)local 
exchange-correlation operators. We label this very simple 
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FIG. 3. (Color online) Experimental and theoretical HOMO- 
LUMO gaps in electronvolts. Red circles: experimental val- 
ues; light blue triangles up: LDA Kohn-Sham HOMO-LUMO 
gap; green squares: non-self-consistent Go Wo (LDA) value; 
black diamonds: GW value with self-consistency on the eigen- 
values; green stars: non-self-consistent GoWo(HFdiag) val- 
ues (see text). The two down-pointing empty triangles are 
B3LYP/6-31G(d) HOMO-LUMO gap values from Refs. [H- 
[55I for Ceo and the H2P porphin. 



scheme GoWo(HFdiag)- This approximation was tested 
by Hahn, Schmidt and Bechstedt^ in the case of three 
small molecules (silane, disilane, water), arguing as we 
do that the Kohn-Sham eigenvalues are too bad a start- 
ing point to evaluate the timc-ordcrcd Green's function 
and the screened potential. Such an approach is also a 
variation on the GoWo(HF) scheme recently introduced 
in Ref. [s^ which was shown to yield the best ionization 
energies for small molecules. With increasing size and 
number of electrons, the part of correlations in the self- 
energy is expected to become more important and using 
Hartree-Fock eigenstates/eigenvalues as a starting point 
for the much larger systems we study may, in principle, 
not be better than using (semi)local functionals for gen- 
erating the starting cigenstates. This is what wc now 
explore. 

For sake of comparison, wc have studied the two 
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HOMO-LUMO gap 




LDA-KS 


GoWo(LDA) 


GW 


GoWo(HFdiag) 


Experiment 


anthracene 


2.25 


6.15 


6.74 


6.86 


6.9" 


tetracene 


1.57 


5.03 


5.58 


5.69 


5.9" 


pentacene 


1.10 


4.21 


4.76 


4.86 


5.2" 


Ceo 


1.58 


4.44 


4.91 


5.08 


4.9" 


MAE 


4.10 (71%) 


0.76 (13%) 


0.22 (3.8%) 


0.10 (2%) 




PTCDA 


1.52 


4.53 


5.0 


5.11 




H2P 


1.94 


4.79 


5.31 


5.44 




H2TPP 


1.82 


4.23 


4.71 


4.91 




HaPc 


1.42 


3.67 


4.03 


4.12 




thiophene 


4.49 


9.93 


10.61 


10.71 




fluorene 


3.59 


7.72 


8.38 


8.54 




benzotliiazole 


3.85 


8.62 


9.40 


9.56 




thiadiazole 


4.29 


10.19 


10.81 


10.89 




benzothiadiazole 


2.94 


7.52 


8.14 


8.23 





TABLE III. HOMO-LUMO gap in eV as obtained from the Kohn-Sham eigenvalues (LDA-KS), non-self-consistent GoWo(LDA) 
calculations, a GW calculation with self-consistency on the eigenvalues (GW), and a non-self-consistent Go Wo calculation 
starting from Hartree-Fock-like eigenvalues {GoWo^HFdiag), see text). MAE is the average mean error in eV for the anthracene, 
tetracene, pentacene and Ceo cases for which experimental band gap data are available. The average error in percent as compared 
to the experiment is indicated in parenthesis. "Ref. |47|. 



small carbon-based conjugated molecules C2H2 and 
C2-ff4 vifhich were investigated by Rostgaard and covifork- 
ers within their full GqWq{HF) scheme. The present 
GoWo(HFdiag) treatment increases the ionization energy 
by 3.48 eV and 3.80 eV for C2H^ and C2H2 respec- 
tively as compared to the LDA values. Such corrections 
compare well with the 3.61 eV and 3.90 eV values ob- 
tained within the full GoWo(-ff-F) scheme of Rostgaard 
and coworkers (as compared to DFT/PBE), emphasizing 
the reliability of the present simplified approximation. 

As compiled in Table II and III (column 5) and in 
Figs. 2 and 3 (green stars), we do find as well that a sin- 
gle shot GoWo(HF(;Mg) calculation provides results which 
are in good agreement with the full GW calculations 
with self-consistency on the eigenvalues. In particular, 
the GoWo(HFdmg) calculations yield much better results 
than the GoWo(LDA) scheme. Such a conclusion agrees 
with that of Rostgaard and coworkers concluding that for 
small isolated molecules, the full GoWo(HF) scheme ac- 
tually outperforms a full self-consistent GW calculation 
where both eigenstates and eigenvalues are updated^. 

Within the present GoWo(IIF£i,;ag) approach, the MAE 
on the ionization energies as compared to experiment is 
0.31 eV, in good agreement with the 0.4 eV result of 
Ref. m for small molecules. Such an agreement indicates 
that the present GoWo(IIFdiag) implementation captures 
most of the features of a full GqWq{HF) approach, sug- 
gesting that LDA and HF eigcnfunctions may not too dif- 
ferent for this set of molecules, a conclusion oftened dis- 
cussed in the literature. Further, the error on the band 
gap, averaged on the calculated values for anthracene, 
tetracene, pentacene and Ceo, for which precise experi- 
mental data are available, is found to be as small as 0.1 
eV (2% error). Such values compare very well with accu- 
rate quantum chemistry calculations with a scheme, the 
GW formalism, which can be applied both to finite size 



and extended systems, and allows to obtain not only the 
band edges, or frontier orbitals, but also the full quasi- 
particle spectrum (see note [sol ) . 

V. CONCLUSIONS 

We have explored the performances of several GW ap- 
proximations for the calculation of the ionization energy 
and HOMO-LUMO gap of thirteen "large" molecules of 
interest for photovoltaic applications, including Cgo, free- 
base porphyrins and phtalocyanine, PTCDA and stan- 
dard donor monomers such as thiophene. Our calcula- 
tions are based on a gaussian-basis implementation with 
full dynamical effects through contour deformation tech- 
niques. Due to the dramatic error on the HOMO-LUMO 
gaps obtained with (semi)local functionals, we find that 
the standard non-selfconsistent CqWo calculations based 
on input LDA eigenstates performs rather poorly, in par- 
ticular in evaluating the HOMO-LUMO gaps. A sim- 
ple self-consistency on the eigenvalues used to build G 
and W provides much better results. As an even sim- 
pler scheme, a non-self-consistent GoWo(HFdiag) start- 
ing from Hartre-Fock like eigenvalues provides equivalent 
results. Both the GW and GoWo(HFdiag) approaches 
provide ionization energies with a mean average error 
within ^ 0.3 eV (^ 4%) of the experiment. Concerning 
the HOMO-LUMO gaps, with a limited number of ex- 
perimental data, the same GW and GoWo(HFdiag) ap- 
proaches yield a mean average error of 0.1-0.2 eV (2-4%), 
in much better agreement than the 4.1 eV (71%) error 
within DFT/LDA, but also in significantly better agree- 
ment than the 0.76 eV (13%) error within the "standard" 
Go Wo (LDA) approach. The possibility of performing 
GW calculations for molecules comprizing several dozens 
of atoms with reasonnable computer time and accuracy. 
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with a scheme allowing to obtain the full quasiparticle 
spectrum of both finite size and extended systems, opens 
the way to the investigation of organic photovoltaic sys- 
tems with techniques that may possibly compete with 
well-established quantum chemical approaches. 
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